# THIS CODES RUNS THE SIMULATIONS FOR THE PAPER
# "Where's the Bias", by Cesar Zucco Jr.
#
# Plots the histograms with the fited beta distributions
# Exemplifies the changes in the beta throughout the range
#
# It simulates the seat shares given the vote share at three levels of dispersion of vote
# shares across districts (sd = approx 0, sd = 10, sd = 20). For each of these variance levels the aggregate 
# vote share is allowed to varry from 10 to 90% with 1% intervals. For each aggregate vote
# share it draws 300 sets of 60 district vote shares from a beta distribution with fixed variance and 
# varying mean. More details on the simulation write up. 
#
#
# Bottomline: As variance in voteshare across districts increases, the system becomes more proportional!! 
# Before running the simulation the code also fits the data to a beta distribution by maximum likelihood

# Fitting the data to a beta distribution
library(foreign)
setwd( "~/Dropbox/Data/Paper-Chilean Electoral System")
chile.data<-read.dta("Temp/chile2.dta")      #Inputs the data set from stata format
attach(chile.data)                      #Makes it a dataframe
house <- subset(chile.data, house == 1) #Filters the data to use only Camara returns
y.cct <- house$ccts/100                 #Gets variable ccts (Concertacion vote share) and calls it y.cct
y.der <- house$ders/100                 #Same for the right

loglikelihood <- function(theta,y) {
   a <- theta[1]
   b <- theta[2]
   llhd <- (a-1)*log(y) + (b-1)*log(1-y) - log(beta(a,b))
   llhd <- sum(llhd)                                        #remember, we're maximizing the sum of the likelihoods
                                                            #alternatively, we can build the sums into the first line
   return(llhd)
}
  
# Optimize the log-likelihood for CONCERTACION
mle.cct <- optim(c(0.1,0.1),         # starting value for theta
    loglikelihood,               # the log-likelihood function
    method="BFGS",               # optimization method
    hessian=TRUE,                # return numerical Hessian
    control=list(fnscale=-1),    # maximize function instead of minimize  
    y=y.cct)                     # the data

cat("MAXIMUM LIKELIHOOD FOR CONCERTACI?N'S BETA DISTRIBUTIOn:\n")
cat("MLE of theta\n")
print(mle.cct$par)
cat("Hessian at MLE\n")
print(mle.cct$hessian)
cat("Standard Error\n")
print(sqrt(diag(solve(-1 * mle.cct$hessian))))

# Optimize the log-likelihood for ALIANZA
mle.der <- optim(c(0.1,0.1),         # starting value for theta
    loglikelihood,               # the log-likelihood function
    method="BFGS",               # optimization method
    hessian=TRUE,                # return numerical Hessian
    control=list(fnscale=-1),    # maximize function instead of minimize  
    y=y.der)                     # the data

cat("MAXIMUM LIKELIHOOD FOR ALIANZA'S BETA DISTRIBUTIOn:\n")
cat("MLE of theta\n")
print(mle.der$par)
cat("Hessian at MLE\n")
print(mle.der$hessian)
cat("Standard Error\n")
print(sqrt(diag(solve(-1 * mle.der$hessian))))

#GRAPHING
par(mfrow=c(1,1),mar=c(3.5,2,1,1), ask=TRUE)
ruler <- seq(0,1,.001)
hist(y.cct                                         #histogram of the data
    , col="light gray"
    , freq = FALSE
    , main = "Concertacion"
    , xlab = ""
    , cex.axis = 1
    , xlim = c(0,1)
    )
   mtext("Vote Share", side = 1, line = 2)
   
lines(ruler,dbeta(ruler,mle.cct$par[1],mle.cct$par[2]),lty=2)    
savePlot(filename="graph8cct", type=c("ps"))
#my.norm <-rnorm(10000,mean=mean(house$ccts/100),sd=sd(house$ccts/100)) Optional, draws a normal instead of a beta
#dnc <- density(my.norm, kernel = c("gaussian")) 
#lines(dnc$x,dnc$y,lty=3)

hist(y.der                                         #histogram of the data
    , col="light gray"
    , freq = FALSE
    , main = "Alianza"
    , xlab = ""
    , cex.axis = 1
    , xlim = c(0,1)
    )
mtext("Vote Share", side = 1, line = 2)
lines(ruler,dbeta(ruler,mle.der$par[1],mle.der$par[2]),lty=2)    
savePlot(filename="graph8der", type=c("ps"))


# THE PROBLEM OF BETAS WITH DIFFERENT MEANS BUT EQUAL VARIANCE
# Funcion transf was obtained from the first and second moments of the beta distribution
# The a' that makes this expression 0 is the a of interest.
# x is the expected value, or average vote share, that will vary in every iteraction
# v is the variance, which is fixed, and obtained from the data (or from teh fitted beta distribution)
# For more details see the write up for the simulation

par(mfrow=c(1,1))
v <- var(y.der)
a <- mle.der$par[1]
b <- mle.der$par[2]
m <- mean(y.der)

for(m in seq(0.1,0.9,0.01)) {
    transf <- function(a,x,v) { 
        -v + (a*(a+1))/((a+((a*(1-x))/x))*(a+((a*(1-x))/x)+1)) - x^2   
    }
    temp <- uniroot(transf, interval=c(0.001,100), tol = 0.0001, x=m,v=v)
    new.a <- temp$root
    new.b <- (new.a*(1-m))/m
    new.v <- (a*b)/((a+b)^2*(a+b+1))
    plot(ruler,dbeta(ruler,new.a,new.b),type="l")
    par(mfrow=c(1,1),mar=c(4,4,1,1), ask=FALSE)
    }

#SIMULATION
par(mfrow=c(3,1),mar=c(3.5,3.5,1,1),ask=FALSE) #mar(bottom, left, top, right)
for(var. in c(1e-05,0.01,0.04)){        #Determines that the procedure will be repeated 
                                        #for these three values of variance
trials <- 300
seat.dist <- matrix(nrow=1,ncol=5)      #creates the matrix that will store these 5 results: m,sd,seats,lower.ci,upper.ci
for(m in seq(0.1,0.9,0.02)) {           #establishes that the average vote share will vary from 0.1 to 0.9
total.seats<-vector(length=0)           #creates a vector to store total number of seats in each trial

    transf <- function(a,x,v) {         #this function was obtained from the first and second moments of the beta, see above
        (a*(a+1))/((a+((a*(1-x))/x))*(a+((a*(1-x))/x)+1)) - x^2 -v          
    }
    temp <- uniroot(transf
        , interval=c(0.01,100000)
        , tol = 0.0001
        , x=m           #x is the mean vote share across districts, that varies in each "for m" loop
        , v=var.)       #v is the standard deviation, that varies in each of the levels set at the "for var." loop                       
    new.a <- temp$root
    new.b <- (new.a*(1-m))/m

for(i in 1:trials){
v.shares <- rbeta(60, new.a,new.b)      #draws 60 vote shares from the a beta distribution with parameters calculated above
s0 <- v.shares < 1/3                    #creates three vectors (60x1)with seats in each district
s1 <- v.shares >= 1/3 & v.shares <= 2/3 #for all districts, one of the three s0,s1,s2 will be 1 
s2 <- v.shares > 2/3                    #and the other two 0, indicating whether 0, 1 or 2 were won
ss <- cbind(s0,s1,s2)                   #concatenates the three vectors into 1 matrix [60x3]
seats <- ss %*% c(0,1,2)                #multiply the matrix by a [3x1]vector 0,1,2. 
                                        #The result is a vector [60x1] of seats won 
total.seats<-rbind(total.seats
    ,sum(seats))                        #adds the total number of seats and appends it to 
}                                       #bottom of vector "total seats", trials times
sort.seats<-sort(total.seats)           #sorts the total number of seats
lower.ci <- sort.seats[(trials*0.025)] #identifies the 5% conf. interval bounds
upper.ci <- sort.seats[(trials*0.975)]
seats.m <- matrix(c(m,round(sqrt(var.),2),mean(total.seats),lower.ci,upper.ci),ncol=5) 
                                        #saves seat share, sd, number of seats, and bounds info
                                        #a [1,5] vector called seats.m
seat.dist <- rbind(seat.dist,seats.m)   #appends seat.m on the bottom of seat.dist
}

plot(seat.dist[,1],seat.dist[,3]/120
    ,pch = 20
    ,ylim = range(0,1)
    ,xlim = range(0,1)
    ,ylab = " "
    ,xlab = " "
    )
    abline(0,1)
    mtext("Average vote share across districts", side = 1, line = 2, cex=0.7)
    mtext("Seat Share", side = 2, line = 2, cex=0.7)
    text(0,0.95, label=paste("St. Deviation = ",round(sqrt(var.),1)), pos=4, lwd=3)
rm(m,s0,s1,s2,ss,seats,total.seats,sort.seats,lower.ci,upper.ci,seats.m)
}

savePlot(filename="variance", type=c("ps"))
